%% 
function [threshold]= FindThreshold(in)
%



v=in(:);
[n,xout]=hist(v,80);
p=n/length(v);
figure(1);plot(xout,p) %for debug

ut=p(1:end)*[1:1:length(p)]';

for k=1:length(p)
    wk=sum(p(1:k));
    uk=(p(1:k)*[1:1:k]')/wk;
    sigma_b(k)=((ut*wk-uk)^2)/(wk*(1-wk));
end

[val,ind]=max(sigma_b);

w0=sum(p(1:ind));
w1=1-w0;
u0=(p(1:ind)*[1:1:ind]')/w0;
u1=(p(ind+1:end)*[ind+1:1:length(p)]')/w1;

threshold=xout(ind);
end